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1.  Introduction 


The  present  and  anticipated  roles  of  tactical  aircraft  impose  serious  challenges  for  the 
design  of  aircraft  exhaust  nozzles.  Present  or  expected  aircraft  requirements  include 
decreased  aftbody  drag  through  improved  airframe-propulsion  integration,  enhanced 
aircraft  maneuver  capability,  short  take-off  and  landing  (STOL),  reduced  aircraft 
observables  and  increased  supersonic  cruise  range  ( [11-  [6]). 

The  accurate  prediction  of  the  performance  of  exhaust  nozzle  flows  is  challenging  due  to 
several  factors,  including: 

•  Strong  viscous-inviscid  interaction:  Flow  separation  in  the  vicinity  of  the 
exhaust  nozzle(s)  is  a  common  characteristic,  due  to  strong  shock-boundary 
layer  interaction  or  adverse  pressure  gradients  in  the  absence  of 
shocksf  [7)-  [10]). 

•  Complicated  geometrical  shapes:  The  integration  of  exhaust  nozzles  into  the 
airframe  results  in  non-simple  three  dimensional  geometrical  shapes  which 
can  generate  complex  3-D  flow  patterns.  For  example,  the  integration  of  twin- 
jet  axisymmetric  nozzles  with  a  rectangular  fuselage  is  sometimes  achieved 
through  the  use  of  boattailed  "gutter"  interfairings,  which  can  adversely  affect 
vehicle  drag  [11]. 

•  Unsteady  flow  fields:  The  geometry  of  certain  exhaust  nozzle-airframe 
integrations  is  possibly  subject  to  low  frequency  unsteady  fluid  motion.  For 
example,  the  region  between  two  widely-separated  twin-jet  axisymmetric 
nozzles  may  be  anologous  to  an  "open"  cavity,  which  are  observed  to  display 
self-sustained  oscillations  in  a  number  of  different  geometrical  configurations 
([12]-  [15]). 

The  current  approach  to  the  design  and  evaluation  of  exhaust  nozzles  relies  heavily  on 
sub-  and  full-scale  model  tests  and  empirical  correlations(  [2],  [4]- [8],  [16]- [18]). 
Theoretical  analysis  typically  consists  of  a  combined  inviscid-boundary  layer  approach, 
with  empirical  corrections  to  account  for  the  discrepancy  between  prediction  and 
experiment  [6],  This  method,  however,  is  incapable  of  handling  strong  viscous-inviscid 
interactions,  and  has  failed  to  accurately  predict  the  performance  of  some  advanced  nozzle 
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configurations  [6].  The  evolution  towards  improved  airframe-propulsion  integration  and  the 
increased  geometric  complexity  in  nozzle  designs  (e.g.,  single  expansion  ramp  nozzle 
[SERN],  2-D  converging-diverging  [C/D]  nozzles  with  yaw  and  pitch  vectoring)  results  in  a 
dramatic  increase  in  the  number  of  configurations  to  be  tested,  and  an  associated  rise  in 
development  costs. 

In  recent  years,  the  capability  for  accurate  and  efficient  numerical  simulation  of  complex 
flows  involving  strong  viscous-inviscid  interactions  has  been  significantly  enhanced  by  two 
factors.  First,  the  advent  of  modern  high-speed  vector-processing  computers  (such  as  the 
CRAY-XMP  and  CYBER  205)  has  afforded  typically  one  to  two  orders  of  magnitude 
improvement  [19,  20,  21, 22]  in  computational  efficiency  compared  to  the  earlier  generation 
computers  such  as  the  CDC  6600  and  IBM  370/168.  Second,  the  development  of  efficient 
implicit  and  hybrid  implicit-explicit  numerical  algorithms  tor  the  full  unsteady  mean 
(Reynolds-averaged)  Navier-Stokes  equations  ( [21],  [23]- [27])  has  also  enhanced 
computational  efficiency.  As  a  consequence,  emerging  design  methodologies  for  exhaust 
nozzles,  as  well  as  other  aircraft  components,  envision  the  utilization  of  full  Navier-Stokes 
numerical  simulations  as  part  of  a  hierarchy  of  theoretical  approaches  which  also  include 
the  traditional  combined  inviscid-boundary  layer  analysis  and  the  recently  developed 
parabolized  Navier-Stokes  methods  [28],  Together,  these  approaches  can  be  combined 
with  experimental  testing  in  order  to  minimize  the  number  of  configurations  in  the  required 
experimental  test  matrix  [29],  thereby  allowing  greater  attention  to  a  selected  number  of 
experimental  test  configurations,  and  a  reduction  in  overall  development  costs. 

A  number  of  steady  2-D/axisymmetric  numerical  simulations  of  nozzle  exhaust  flow 
fields  have  been  performed  in  recent  years  using  the  compressible  Navier-Stokes 
equations.  These  include  the  computations  of  axisymmetric  nozzles  by  Mikhail,  Hankey 
and  Shang  [30]  and  Hasen  [31]  and  two-dimensional  nozzles  by  Cline  and  Wilmoth  [32] 
and  Perry  [33]  Steady  three-dimensional  exhaust  nozzle  computations  have  been 
performed  using  the  mean  Navier-Stokes  equations  However,  unsteady  two-dimensional 
compressible  laminar  Navier-Stokes  calculations  have  been  performed  lor  spike-tipped 


bodies  by  Shang,  Smith  and  Hankey  [15]  and  for  flow  past  a  cylinder  by  Shang[34],  in 
addition,  unsteady  2-D  turbulent  Navier-Stokes  simulations  have  been  computed  by 
Levy  (35)  for  a  circular  arc  transonic  airfoil,  and  by  Steger  and  Bailey  [20]  for  the  transonic 
flow  past  the  F-80  wing.  In  these  latter  cases,  the  unsteady  flow  was  characterized  by  low 
frequency  motion  (time  scales  on  the  order  of  the  mean  flow),  and  the  effects  of  the  high 
frequency  turbulent  motion  was  incorporated  through  an  algebraic  turbulent  eddy  viscosity 
model. 


The  primary  objective  of  this  research  is  to  develop  and  evaluate  the  ability  to  simulate 
complex  nozzle  flows.  The  solution  of  problems  with  finite  difference  methods  (CFD) 
consists  of  three  major  phases:  mesh  generation,  flow  computation  and  flow  analysis.  The 
complex  geometries  encountered  in  nozzle  flows  consist  of  curved  or  irregularly  shaped 
boundaries  for  which  grid  generation  becomes  a  task  in  itself.  Uniformly  discretized  grids  in 
the  physical  domain  are  inadequate  since  in  addition  to  difficulties  encountered  in 
application  of  boundary  conditions,  the  accuracy  of  the  computation  may  be  affected  and  it 
may  be  impossible  to  achieve  adequate  resolution  in  regions  with  large  gradients  with  a 
limited  number  of  mesh  points.  A  number  of  techniques  of  grid  generation  for  finite 
difference  applications  are  summarized  in  [36]  and  [37],  From  a  survey  of  the  literature, 
though  a  number  of  techniques  have  been  proposed  for  3-D  grid  generation,  evidently  the 
focus  of  application  has  been  on  2-D  cases.  A  need  further  exists  to  incorporate 
interactive  graphics  [38]  into  the  grid  generation  process  to  reduce  the  time  required  to 
generate  complex  grids.  This  report  describes  the  implementation  of  one  particular 
method  (the  multisurface  technique  of  Eiseman  [39]). 

The  focus  of  the  nozzle  flow  simulation  is  on  steady  flow  in  a  nonaxisymmetric  wedge 
nozzle  at  a  freestream  Mach  number  of  1 .2.  The  choice  of  this  case  is  dictated  by  the 
existence  of  experimental  investigations  with  surface  pressure  measurements  (subsection 
3.2).  The  explicit-implicit  algorithm  employed  (subsections  3.4  and  3.5)  is  applicable  for 
axisymmetric  nozzles  as  well,  although  some  additional  work  may  be  necessary,  however, 
to  resolve  coordinate  transformation  singularities. 
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The  application  of  the  multi-surface  technique  of  Eiseman  is  described  in  Section  2. 
Section  3  describes  the  nozzle  flow  computation  including  the  theoretical  model  for  nozzle 
flow  (governing  equations  and  turbulence  model),  grid  details  and  preliminary  numerical 
results.  Section  4  concludes  with  some  remarks  on  the  current  ability  to  simulate  nozzle 
flows  and  identifies  areas  of  future  work. 

2.  Three-Dimensional  interactive  Grid  Generation  Using  the 
Multisurface  Technique 

2.1.  Introduction 

The  use  of  body-fitted  grids  considerably  simplifies  the  problem  of  flow  simulation  using 
finite  difference  techniques  in  domains  with  curved  boundaries.  For  such  domains,  the 
employment  of  Cartesian  and  other  standard  meshes  becomes  cumbersome  since 
interpolation  is  generally  necessary  at  the  boundaries.  It  is  also  difficult  if  not  impossible  to 
concentrate  the  grid  in  arbitrary  regions  of  space  without  introducing  a  large  number  of 
unnecessary  mesh  points  [40].  The  basic  concept  is  to  map  the  given  complicated  physical 
domain  and  the  governing  equations  into  a  topologically  simpler  domain  as  shown  in  Fig. 
2-1,  which  demonstrates  the  procedure  for  flow  about  an  airfoil.  A  significant  amount  of 
effort  has  been  focussed  recently  on  the  development  of  adequate  grid  generation 
techniques  applicable  to  general  geometries.  Ref.  [41]  provides  an  overview  of  this  effort. 

The  standard  approach  to  grid  generation  involves  the  application  of  a  selected  grid 
generation  algorithm  (e.g.,  the  multisurface  technique  of  Eiseman  [39])  to  a  specific 
configuration.  The  algorithm  may  utilize  one  or  more  techniques  for  grid  point  distribution 
(e  g.,  the  location  of  intermediate  surfaces  in  the  algebraic  grid  generation  method  of 
Eiseman  [39]).  The  grid  is  subsequently  examined  for  satisfaction  of  the  specified  grid 
requirements  (e  g.,  orthogonality  of  the  grid  lines  near  solid  boundaries).  If  the  grid  is 
deemed  unsatisfactory,  the  input  to  the  grid  generation  algorithm  is  modified  (e  g.,  the 
location  of  the  intermediate  surfaces),  and  the  grid  recomputed.  This  process  of  grid 
generation,  diagnosis,  and  modification  of  input  variables  is  continued  until  a  satisfactory 


4 


T)  =  T|  (r,y) 

(*.y) 


x  =  x  (n,S ) 
y  =  y  (n4  ) 


l 


Figure  2-1 :  Curvilinear  Grids  -  Physical  and  Transformed  Planes 
grid  is  achieved. 

The  specific  grid  requirements  of  a  given  application,  however,  may  not  be  completely  or 
easily  met  by  the  above  approach.  In  the  current  research  effort,  a  two-phase  methodology 
is  developed  (Fig.  2-2).  Phase  I  employs  a  grid  generation  algorithm  (the  multisurface 
technique  of  Eiseman  [39])  and  an  algorithm  for  control  of  grid  point  distribution  (the 
method  of  Smith  et  at [42]).  The  purpose  of  the  first  phase  (Phase  I)  is  to  obtain  an 
approximate  grid  satisfying  part  of  the  grid  requirements.  Phase  II  employs  graphical 
techniques  to  directly  modify  local  areas  of  the  grid  to  meet  the  remaining  grid 
requirements.  The  second  phase  is  performed  independently  of,  and  subsequent  to,  the 
grid  generation  by  the  specified  algorithm  (Phase  I),  Both  phases  are  incorporated  in  an 
interactive  environment  utilizing  color  graphics. 

It  is  evident  from  a  survey  of  the  literature  that  while  significant  effort  has  been  focussed 
on  the  generation  of  2-D  grids,  there  have  been  relatively  fewer  applications  of  3-D 


5 


Figure  2-2:  Grid  Generation  Approach  of  this  Research 
methods.  This  is  because  the  volume  and  visual  complexity  of  data  specification  increases 
significantly  making  the  use  of  interactive  graphics  (in  particular  analog  graphical  input) 
imperative.  While  the  importance  of  graphics  for  evaluation  of  generated  grids  has  long 
been  recognized,  the  applicability  of  interactive  graphics  techniques  in  the  grid  generation 
phase  itself  has  not  been  extensively  explored  (with  the  notable  exception  of  the  work  of 
Smith  ef  a/ [42]). 

This  section  describes  a  three-dimensional  graphics-based  interactive  grid  generation 
program.  The  use  of  graphics  is  extended  to  ail  stages  of  grid  generation  with  particular 
emphasis  on  direct  modification.  The  code  has  the  capability  of  generating  meshes  for 
general  configurations.  Two  typical  examples  are  presented 
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2.2.  Grid  Requirements 

While  the  intent  of  this  research  is  to  develop  a  method  with  broad  application,  for 
concreteness  the  twin-jet  axisymmetric  nozzle  is  considered  (see  Fig.  2-3).  This 
configuration  displays  several  of  the  principal  features  related  to  grid  topologies  of  practical 


interest.  These  are  summarized  as  follows: 

•  Body  shape  exhibits  transition  in  shape 

•  Two  flows,  inner  and  outer,  which  meet  at  the  nozzle  exit  exist.  The  mixing  of 


these  constitutes  the  wake  region  which  must  be  resolved  with  an  adequate 
grid. 

In  addition  to  these,  the  numerical  algorithm  employed  in  the  solution  of  the  flow 


& 
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equations  introduces  some  constraints.  We  arrived  at  the  following  general  constraints  by 
considering  a  number  of  algorithms  -  the  explicit  method  of  McCormack  [43],  the  hybrid 
explicit-implicit  method  of  Knight  [44]  and  the  implicit  method  of  Beam  and  Warming  [25], 

•  Orthogonality  is  desirable  in  boundary  layer  regions,  especially  for 
implementation  of  algebraic  eddy  viscosity  models. 

•  Precise  point  distribution  control  is  necessary  to  resolve  viscous  regions  of  the 
boundary  layer  and  the  inviscid  regions  outside  the  boundary  layer. 

To  simplify  the  graphical  input  requirements,  the  grids  are  generated  in  axis-normal 
sections  and  stacked  together.  The  control  method  described  below  is  employed  to 
transmit  point  distribution  information  in  the  axial  direction  to  prevent  excessive  skewing  of 
grid  lines. 

2.3.  Phase  I:  Grid  Generation  by  Multisurface  Interpolation 

The  Phase  I  grid  is  generated  by  applying  the  multisurface  technique  algorithm 
described  by  Eiseman  in  [39],  In  this  method,  a  number  of  parameterized  surfaces  are 
introduced  as  shown  in  Fig.  2-4.  The  basic  formula  [39]  is  : 

AM 
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Figure  2-4:  The  Muitisurface  Technique 


where 


Gt(r)= j'  vt(x)dx. 
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P,  through  PN  are  the  parameterized  surfaces  (with  common  parameter  ?),  r  is  the 
spanning  variable  (takes  the  values  r,  and  rN  ,  at  the  inner  and  outer  bounding  surfaces 
respectively),  the  tunctions  y  are  interpolants.  Choosing  Aw  as  in  equation  3  leads  to  the 
bounding  surfaces  being  part  of  the  transformation  Depending  upon  the  choice  of  function 
for  yK,  two  types  of  interpolation  -  global  and  local  may  be  specified  employing  respectively 
globally  defined  polynomials  in  r  and  piecewise  linear  functions  non-vanishing  only  in  a 
local  region.  Both  methods  are  incorporated  in  the  current  research. 

Since  precise  grid  point  control  through  iterative  specification  of  the  intermediate 
surfaces  tends  to  be  tedious  even  in  an  interactive  environment,  the  control  method  of 
Smith  et  al  [42]  is  employed.  This  control  method  provides  non-uniform  discretization  of 
any  general  physical  variable  in  a  convenient  interactive  graphical  environment.  A  few 
points  are  first  digitized  on  a  plot  whose  ordinate  is  a  parametric  variable  trom  the  physical 
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plane  and  whose  abscissa  is  a  normalized  variable  from  the  computational  plane  A 
smooth  cubic  spline  is  then  passed  through  these  points  Shallower  regions  along  this 
curve  correspond  to  higher  concentration  of  the  physical  variable  discretization 

An  example  of  the  application  of  the  control  method  in  this  research  is  shown  in  Fig  2-5. 
The  first  step  in  the  control  of  the  point  distribution  in  Fig.  2-5(a)  is  to  plot  the  cumulative 
arc  length  profile  (CALF)  along  a  cubic  spline  through  the  set  of  points  against  point 
number  (curve  marked  "initial"  in  Fig  2  5(b))  A  new  CALP  (marked  "final")  is  then 
digitized  or  specified  analytically.  In  the  present  case,  the  exponential  formula 


c^-  1 

r  =: - -  ;  0  <  r  <  I 

ek- I 


(4) 


is  employed  [39]  with  a  k  value  of  1  0.  The  new  CALP  may  also  be  specified  to  be  similar 
to  CALPs  of  other  distributions.  The  cubic  spline  coefficients  associated  with  the  set  of 
points  in  Fig.  2-5(a)  are  then  used  in  conjunction  with  the  new  CALP  to  obtain  the  new 
point  distribution  shown  in  Fig.  2-5(c).  Note  that  since  the  "final"  curve  in  Fig  2-5(b)  has  a 
smaller  value  of  slope  near  point  1,  the  final  distribution  has  a  correspondingly  higher 
concentration  near  point  1 . 

2.4.  Phase  II:  Direct  Grid  Modification 

This  phase  incorporates  a  series  of  powerful  graphics  based  options  that  facilitate  grid 
modification  influencing  either  the  entire  grid  or  only  points  along  a  specific  grid  line  This 
phase  is  particularly  useful  in  obtaining  special  local  grid  characteristics  (e  g  ,  corners) 
without  adversely  affecting  the  entire  grid  (see  subsection  2.6)  For  example,  since  the 
multisurlace  technique  guarantees  only  that  the  bounding  surfaces  are  part  of  the 
transformation,  it  is  necessary  to  employ  special  techniques  to  obtain  precise  transverse 
boundary  shapes  (e.g.,  the  extreme  g  -  constant  lines  in  Fig  2-4  in  the  physical  domain) 

The  application  of  Phase  II  generally  begins  with  the  identification  of  a  set  of  points 
along  a  grid  line.  The  points  may  be  relocated  in  a  variety  of  ways  Some  important  options 
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•  Rotation:  The  points  may  be  rotated  about  any  axis  perpendicular  to  the  plane 
of  the  grid.  The  point  ot  intersection  between  the  axis  and  the  plane  may  be 
specified  either  by  digitizer  or  by  keyboard  input  When  the  entire  grid  is 
identified,  as  opposed  to  a  set  ot  points  along  a  grid  line,  rotation  may  be 
affected  about  an  arbitrary  axis  in  three  dimensional  space.  This  permits  rapid 
generation  of  grids  for  solids  of  revolution  Reflection  about  any  axis  may  be 
achieved  by  a  180°  rotation  about  the  appropriate  axis 

•  Horizontal  or  vertical  alignment  All  identified  points  may  be  aligned  either 
vertically  or  horizontally.  The  use  of  this  option  is  illustrated  later  for  the  twmjet 
nozzle  case, 

•  Modification  with  cubic  splines  1  tie  identified  points  may  be  redistributed  with 
the  process  described  earlier  m  Fig  2-5 

•  Point  position  modification  Fach  point  may  be  relocated  with  graphical  or 
alphanumeric  input. 

•  Change  in  total  number  ol  points  The  total  number  of  points  in  either  direction 
may  be  increased  or  reduced  without  affecting  the  relative  concentration 
characteristics  This  is  done  by  first  computing  the  cubic  spline  coefficients 
and  CALP  for  each  grid  line  and  inverting  the  CALP  at  the  new  number  of  grid 
points  to  obtain  the  new  grid  point  distribution 

2.5.  Description  of  Code 

The  three  dimensional  graphics-based  code  is  written  in  Fortran  77  and  utilizes  graphics 
device  driver  routines  specific  to  the  Tektronix  41 XX  series  from  the  Tektronix  Interactive 
Graphics  Library  (IGL).  The  program  is  implemented  on  a  Prime  9955  computer.  Input  to 
the  code  is  provided  interactively  in  either  a  digital  (keyboard)  or  analog  (graphics  cursor) 
fashion  Analog  input  permits  the  specification  of  points  in  physical  coordinates  such  as 
necessary  to  define  intermediate  surfaces  and  CALPs.  identify  grid  points  for  direct 
processing  and  window  areas  for  visual  diagnostics  The  sequence  of  execution  of  the 
various  modules  may  be  determined  by  the  user  Tins  is  called  the  "command  driven" 
execution  mode  and  provides  significant  flexibility  in  the  grid  generation  process  Complete 


details  of  the  code  are  provided  in  [45] 

The  grid  generation  code  is  split  into  five  modules  These  are  :  surface  definition,  phase 

1  grid  generation,  phase  II  grid  generation,  diagnostics  and  grid  structuring 

(a)  Surface  definition:  This  is  generally  the  first  step  of  the  grid  generation  process  and 
results  in  the  definition  on  a  set  of  surfaces  (inner,  outer  and  intermediate)  The  surfaces 
may  be  defined  interactively  in  one  of  a  number  of  ways  such  as  through  the  "joy  disk" 
digitizer,  disk  files,  use  of  standard  ("superellipse")  profiles,  normals  to  any  previously 
defined  surfaces  or  on  rays  passing  through  corresponding  points  of  a  pair  of  surfaces  It  is 
also  possible  to  introduce  point  singularities.  Once  the  shape  of  the  surface  is  defined,  the 
point  distribution  on  these  surfaces  may  be  adjusted  with  the  control  method  described 
earlier. 

(b)  Phase  I  grid  generation:  In  this  module,  the  formulas  of  subsection  2  3  are  employed 
to  generate  an  initial  grid  The  spanning  variable  r  may  be  discretized  in  a  non-uniform 
fashion  using  the  control  plane  This  is  not  effective  means  of  control  In  the  applications  to 
be  described,  a  uniform  discretization  is  employed 

(c)  Phase  II  grid  modification  This  module  implements  the  options  outlined  in  subsection 

2  4. 

(d)  Diagnostics:  A  number  of  interactive  diagnostic  measures  are  incorporated  These 
include  window  options  which  permit  detailed  telescopic  examination  of  the  grid  Two 
quantities  indicative  of  grid  structure,  the  Jacobian  and  Skewness  measure  may  also  be 
monitored.  The  skewness  measure  (SM)  is  defined  as 

l  k 

SM  =  £  |W'  <),|  (5) 

k  ii 

where  k  is  the  number  of  angles  surrounding  a  point  and  the  n,  are  the  angles  measured  in 
degrees 

(e)  Grid  structuring:  This  is  generally  the  last  executed  module  At  any  given  section. 
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grids  may  be  generated  in  more  than  one  zone  These  zones  are  combined  to  form  a 
composite  grid  through  the  grid  structuring  process 


The  code  generates  grids  meeting  all  the  requirements  outlined  in  subsection  2  2  The 
physical  features  of  most  geometries  of  interest  are  adequately  treated  with  the  above 
surface  definition  techniques.  Orthogonality  at  the  boundary  is  achieved  by  defining  the 
first  intermediate  surface  along  normals  to  the  inner  boundary  Precise  point  distributions 
are  obtained  with  the  control  method  of  subsection  2  3 

The  general  strategy  is  to  define  intermediate  surfaces  to  satisfy  £  grid  line  shape  and 
point  distribution  and  q  line  shape  alone  to  obtain  an  initial  grid  (see  Fig  2-4).  The  graphics 
based  a  posteriori  modification  technique  (Phase  II)  is  then  employed  to  obtain  the  precise 
required  distribution  along  the  q  lines  The  locations  of  the  numerous  stacked  sections 
control  grid  concentration  in  the  remaining  third  direction 

2.6.  Applications 

Two  applications,  chosen  to  display  the  versatility  of  the  method,  are  illustrated.  The  first 
case  is  the  twinjet  axisymmetric  nozzle  (Fig  2  3).  while  the  second  represents  an  aircraft 
section  with  an  underbelly  inlet  In  both  applications,  though  grids  were  generated  with 
both  global  and  local  interpolation  formulas,  only  those  generated  with  local  interpolation 
are  discussed  since  they  possess  superior  boundary  orthogonality  characteristics  Only  a 
few  representative  grids  are  shown  with  the  understanding  that  the  characteristics  may  be 
altered  using  the  methods  outlined  earlier.  Details  of  the  diagnostic  phase  are  excluded  for 
the  purpose  of  brevity. 

(a)  Twin-jet  axisymmetric  nozzle  The  general  case  is  a  twin-jet  nozzle  with  an 
interfairing  whose  trailing  edge  is  of  finite  thickness  and  located  upstream  of  the  nozzle  exit 
plane.  The  general  grid  configuration  at  an  axial  station  upstream  of  the  nozzle  is  shown  in 
Fig  2-6  An  appropriate  topological  plane  is  shown  in  Fig  2-7  Since  the  shaded  areas 
represent  solid  regions,  no  flow  computations  need  be  performed  here  Note  that  the  entire 
line  AOUNM  maps  on  to  the  center  of  the  nozzle  and  as  such  this  represents  a  singularity 


Upstream  of  Interfairing  Trailing  Edge 

Transformed  Plane 


4 


Figure  2-7:  Transformed  Plane  for  Twinjet  Nozzle 


Figure  2-8:  Sub-Zones  at  Typical  T  winjet  Nozzle  Section 
has  a  vertical  segment  (EF)  to  facilitate  application  of  symmetry  boundary  conditions  there 
If  the  surface  EFG  is  directly  employed  in  the  multisurface  technique,  the  slope 
discontinuity  propagates  into  the  field  This  is  an  undesirable  situation  since  it  may 
adversely  affect  the  flow  computation  In  this  case,  the  multisurlace  technique  is  actually 
employed  with  a  fictitious  elliptic  outer  surface  (see  Fig  2-8)  which  is  later  stretched  out 
This  is  a  typical  example  where  the  application  of  the  horizontal  or  v<  dical  alignment  option 
(Phase  II)  is  useful  (see  subsection  2  4)  Zone  3  is  a  special  zone  since  the  grid  point 
distributions  along  two  adjacent  sides  (QC  and  EQ)  are  known  from  Zones  1  and  2 
respectively.  The  third  boundary  ED  must  be  vertical  to  conform  to  the  vertical  line  of 
symmetry.  The  ordinates  of  the  points  along  this  line  are  chosen  to  be  identical  to  those  for 
curve  QC  and  the  grid  is  generated  by  simple  algebraic  linear  interpolation  There  are  no 
intermediate  surfaces  for  this  zone  and  the  multisurface  technique  is  never  applied 

The  actual  intermediate  surface  distributions  for  Zones  i  and  2  are  shown  in  F  ig  2-10 
Fig.  2-11  displays  the  subgrids  and  Fig  2-12  shows  the  final  composite  grid  obtained  by 
structuring  the  grid  sections  of  Fig  2  11  All  grids  are  exponentially  smoothed  towards 
their  respective  inner  boundaries  (see  equation  <4)  Other  distributions  may  be  rapidly 
generated  with  the  precise  control  method  described  earlier  Grids  for  other  sections 


Zone  1 


Inner  boundary:  Singularity  corresponding 
to  center  of  nozzle  (A). 

Outer  boundary:  Nozzle  Surface  (CQV). 
Intermediate  Surfaces  (1):  Circular  (XY). 
Points:  20  in  eta  direction. 

15  in  zeta  direction. 


Inner  boundary:  Nozzle  surface  (QV). 

Outer  boundary:  Fictitious  ellipse  (EG). 
Intermediate  Surfaces  (2):  Along  normals 

and  midway  between  inner  and  outer 
respectively. 

Points:  15  in  eta  direction 

20  in  zeta  direction 


Three  side  distributions  known: 
EQ  -  from  Zone  2. 

QC  -  from  Zone  1 . 

ED  -  from  QC. 

No  intermediate  surfaces. 


Zone  3 


Figure  2-9:  Intermediate  and  Bounding  Surface  Description  tor 

Twiniet  Nozzle 


upstream  and  downstream  of  the  nozzle  exit  may  be  generated  in  a  similar  fashion 


(b)  Aircraft  section  with  underbelly  inlet:  The  grid  for  this  configuration  is  generated  in 
five  zones  shown  in  Fig.  2-13  The  proportions  employed  are  chosen  to  retain  the  relevant 
generic  features  of  such  a  cross  section.  Zones  1  through  4  represent  regions  of  external 
flow.  Zones  2  and  3  are  reflections  of  Zones  1  and  4  respectively  about  a  vertical  axis 
Zone  5  corresponds  to  the  inlet  region  at  any  section.  Fig.  2-14  displays  the  intermediate 
sudace  set  for  each  zone.  The  transverse  boundaries  ABC  and  FBC  are  made  to  confirm 
precisely  to  the  solid  boundary  with  Phase  II  techniques  discussed  in  subsection  2  4  Fig 
2-15  shows  the  exponentially  smoothed  final  grid. 

2.7.  Conclusion 

A  3-D  grid  generation  code,  based  on  the  multi  surface  technique  of  Eiseman  and  the 
control  method  of  Smith  et  al  has  been  developed  The  code  is  executed  interactively  and 
strongly  emphasizes  the  use  of  color  graphics  User  input  may  be  both  digital  (i.e  .  from  a 
file  or  terminal  keyboard)  and  analog  (i.e  ,  a  graphics  cursor  or  "joydisk").  The  3-D  grid 
system,  which  is  generated  by  "stacking"  a  series  of  2-D  grid  systems,  may  be  viewed 
graphically  with  provision  tor  enlargement  ("zoom")  of  specific  sections  of  interest. 

The  grid  system  may  be  modified  using  a  variety  of  techniques  which  permit 
redistribution  of  grid  points  and  redefinition  of  boundary  surfaces.  The  modification  is  fully 
graphics  based  and  provides  a  variety  of  grid  diagnostics  including  graphical 
representation  of  the  Jacobian,  angle  of  intersection  of  grid  lines  and  distribution  of  grid 
points  versus  distance  along  any  specked  coordinate  line 

The  various  features  of  the  code  may  be  "command  driven"  by  the  user.  Consequently, 
the  user  may  proceed  to  develop,  diagnose  and  modify  the  grid  system  as  desired 
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Zone  1 


Inner  boundary:  Circular  (AD). 

Outer  boundary:  Circular(EC). 
Intermediate  surfaces  (4): 

Points:  20  in  circumferential  direction. 
30  in  radial  direction. 

Zone  3:  Reflection  of  Zone  1  about  DE. 


Inner  boundary:  Inlet  exterior  (FG). 
Outer  boundary:  Circular  (CH). 
Intermediate  surfaces  (4): 

Points:  20  in  circumferential  direction. 
30  in  radial  direction. 

Zone  4:  Reflection  of  Zone  2  about  GH. 


Zone  2 


Zone  5 


Inner  boundary:  Singularity  corresponding 
to  center  of  inlet  (I). 

Outer  boundary:  Elliptic  (JJ’J). 
Intermediate  Surfaces  (1):  Elliptic 
Points:  17  in  circumferential  direction. 

10  in  radial  direction. 


Figure  2-14:  Intermediate  Surface  Structure  tor  Aircraft  Section 
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Figure  2-15:  Final  Grid  for  Aircraft  Section 


3.  Simulation  of  Nozzle  Flows 


3.1.  Introduction 

The  principal  objectives  of  this  research  are  to  develop  and  evaluate  the  ability  to 
simulate  complex  nozzle  flows.  One  of  the  fundamental  characteristics  of  such  flows  is  the 
existence  of  two  generally  parallel  streams  of  flows  -  internal  and  external  -  which  meet 
downstream  of  some  trailing  edge  to  form  a  mixing  layer,  the  properties  of  which  are  of 
some  significance.  Of  the  two  generic  types  of  nozzle  configurations  -  axisymmetric  and 
non-axisymmetric  -  the  latter  category  is  chosen  as  the  focus  of  this  research.  This  choice 
is  governed  by  the  existence  of  experimental  data  base  for  verification  purposes.  We 
emphasize  that  the  complexity  of  the  resultant  flow  structure  for  both  axisymmetric  and 
non-axisymmetric  nozzles  is  expected  to  be  similar. 

In  subsection  3.2,  a  brief  description  of  the  experiments  of  Mason  and  Abeyounnis  [46] 
is  provided.  Fig.  3-1  is  a  picture  of  the  configuration  under  consideration.  Since  the 
experiments  employed  for  validation  display  vertical  and  horizontal  symmetry  (zero  angle 
of  attack  and  yaw),  this  assumption  is  integral  to  the  computations.  The  flow  is  computed  in 
one  quadrant  of  the  nozzle  only  (region  ABCDEFGH  in  Fig.  3-1).  A  topological  discussion 
follows  in  subsection  3.3  which  describes  the  physical  to  computational  plane  mapping. 
The  full  Reynolds  averaged  Navier-Stokes  equations  in  mass  averaged  variables  and  two 
different  turbulent  eddy  viscosity  models  for  different  regions  of  the  flow  constitute  the 
theoretical  model  and  are  described  in  subsection  3.4.  The  discretization  and  numerical 
solution  of  the  governing  equations  is  outlined  in  subsection  3.5.  Some  details  of  the 
implementation  of  the  eddy  viscosity  models  in  the  flow  field  are  provided  in  subsection 
3.6.  Subsection  3.7  describes  the  coordinate  transformation  employed  (grid  structure)  and 
certain  geometrical  simplifications  made  so  that  the  problem  is  more  amenable  to  solution. 
A  brief  description  of  the  boundary  conditions  employed  follows  in  subsection  3.8. 
Numerical  details  such  as  the  CPU  time  required  and  history  of  computation  are  discussed 
in  subsection  3.9.  In  subsection  3.10,  the  numerical  results  are  compared  with 
experimental  measurements.  These  comparisons  are  based  on  static  pressure 
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measurements.  A  preliminary  analysis  of  the  flow  field  is  discussed  in  subsection  3.11. 

3.2.  Brief  Description  of  Experiments 

The  experiments  employed  for  comparison  determine  the  geometry  and  flow  conditions 
under  consideration.  These  experiments  were  performed  by  Mason  and  Abeyounis  [46]  in 
the  Langley  16-Foot  Transonic  Tunnel  (Fig.  3-2)-  a  single-return,  continuous-flow, 
atmospheric  wind  tunnel  with  an  octagonal  test  section.  The  experiments  are  performed  at 
wind-off  conditions  for  a  range  of  freestream  Mach  numbers  from  0.6  to  1.2  on  two  non- 
axisymmetric  wedge  nozzles.  For  each  Mach  number,  two  expansion  ratios,  high  (1.20) 
and  low  (1.06),  are  documented.  The  focus  of  this  research  concerns  the  only  supersonic 
case  documented  (Mach  number  =  1.2)  in  the  low  expansion  configuration  described.  The 
high-expansion-ratio  configuration  was  not  investigated.  The  flow  conditions  are  provided 
in  Table  3-1.  The  nozzle  consists  of  five  components:  two  sidewalls,  an  upper  and  a  lower 
flap  and  a  wedge  centerbody.  Side  and  top  views  of  the  internal  and  external  geometry  of 
the  nozzle  configuration  are  provided  in  Fig.  3-3.  This  figure  also  provides  the  precise 
dimensions  employed  in  the  numerical  simulation.  Certain  geometrical  modifications, 
which  are  also  marked,  (subsection  3.7)  are  made  so  that  the  problem  is  more  amenable 
to  solution.  The  nozzle  configuration  is  instrumented  with  14  rows  of  static-pressure 
orifices,  the  locations  of  which  are  described  in  subsection  3.10.  Further  details  of  the 
model  instrumentation  and  techniques  used  may  be  found  in  [46]  and  will  be  omitted  here 
for  the  purposes  of  brevity 


3.3.  Topological  Discussion 

This  subsection  establishes  certain  conventions  and  definitions  which  help  in  the 
description  of  the  implementation  of  the  theoretical  model  and  numerical  algorithm.  For  the 
purposes  of  the  subsequent  discussion,  the  general  regions  of  the  flow  field  are  shown  in 
Fig.  3-4  which  displays  the  top  and  side  views  (the  streamwise  direction  is  from  left  to  right 
in  both)  of  the  nozzle.  Five  boundary  layers  may  be  defined.  Boundary  layers  1  through  3 
constitute  parts  of  the  internal  flow  and  exist  on  the  flap  inner  wall,  inner  sidewall  and  on 
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External  Flow: 


Ttt(K)  P,(kPa)  NRe(/m) 


1.20  339.03  101.04  12.92x106 


Tw  =  331 .4  K  =  Adiabatic  Wall  Temperature 

Internal  Flow: 


NPR  =  2.0 

T,  j  =  300.2  K 

Upstream  Mach  No.  =  0.19 

P, i  =  83.3  kPa 

Legend: 

M  -  Freestream  Mach  number 

Ttt  -  Average  tunnel  total  temperature 

Pt  -  Average  tunnel  total  pressure 

T,  j  -  Average  jet  total  temperature 

NRe  -  Reynolds  number 

Tw  -  Wall  Temperature 

NPR  -  Nozzle  Pressure  Ratio  (Ptj/P  ) 

'  oo 

the  wedge  respectively.  Boundary  layers  identified  4  and  5  occur  in  the  external  flow  on  the 
flap  outer  wall  and  the  outer  sidewall  respectively.  The  mixing  layer  is  formed  when  the 
internal  and  external  flows  meet  at  the  trailing  edge  of  the  flap.  Under  the  assumptions  of 
subsection  3.7,  there  is  no  mixing  layer  downstream  of  the  sidewall. 

The  employment  of  body-fitted  grids  and  genera!  coordinates  permits  the  entire  nozzle 
(including  internal,  external  and  mixing  layer  flows)  to  be  mapped  onto  a  unit  cube.  This  is 
depicted  schematically  in  Fig.  3-5.  The  %  (index  i)  direction  corresponds  to  the  generally 
streamwise  (x)  direction,  the  t)  coordinate  increases  in  the  generally  spanwise  (z)  direction 
and  the  C  coordinate  is  arranged  so  that  the  coordinate  system  2;-tK  is  right-handea  with  y 
increasing  in  the  general  C  direction.  The  index  k  corresponding  to  the  interior  flap  surface 
is  denoted  knozzl8,  while  that  for  the  shell  is  denoted  kshe)(.  Similarly,  the  index  j 
corresponding  to  the  interior  sidewall  surface  is  denoted  j  ,  and  jsheN  is  employed  for 
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(b)  Flap  geometry. 


Figure  3-3  (contd.):  Internal  and  External  Geometry.  From  [46] 
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a)  Top  View 


B.L.  *  4 


Figure  3-4:  General  Flow  Structure  -  Top  and  Side  Views 
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Ccnteriine  of  N'ot/Je 


b)  Computational  Domain 


Figure  3-5:  Physical  and  Transformed  Domains  -  Nonaxisymmetric  Nozzle 


the  exterior  sidewall  surface.  By  construction  (subsection  3.7), 

j shell  ~  Jno2zle+^ 

(6) 

kshell  ~  *nozz/<+1 

In  Fig.  3-5  the  upstream  boundary  is  ABCD,  a  part  of  which  (height  AM,  width  MN) 
corresponds  to  internal  flow.  The  downstream  boundary  is  EFGH  and  is  located  some 
distance  downstream  of  the  wedge  trailing  edge.  Details  of  dimensions  are  provided  in  a 
later  subsection.  The  plane  of  horizontal  symmetry  is  denoted  ADHE,  a  part  of  which 
(UVWX)  corresponds  to  the  upper  surface  of  the  wedge.  ABFE  is  the  vertical  plane  of 
symmetry.  The  freestream  boundaries  (BFGC  and  CDHG)  are  located  far  enough  away 
from  the  nozzle  for  the  application  of  freestream  boundary  conditions.  The  flap  outer  and 
inner  surfaces  are  mapped  on  the  planes  IJKL  and  MNOP  respectively.  For  purposes  of 
clarity,  only  two  adjacent  boundaries  (QR  and  RS)  of  the  external  sidewall  are  shown  in  the 
figure.  The  sidewall  extends  beyond  the  flap  trailing  edge  up  to  the  wedge  trailing  edge  as 
is  clear  from  the  figure  of  the  nozzle  (Fig.  3-3).  Some  assumptions  on  the  flap  trailing  edge 
and  the  sidewall  are  described  later. 


The  body-fitted  grid  is  generated  in  axis-normal  (E,  =  constant)  planes  for  simplicity.  A 
general  physical  plane  and  the  corresponding  computational  plane  is  shown  in  Fig.  3-6.  At 

any  cross  section,  the  following  naming  convention  is  employed. 

•  Region  1  :  internal  flow  (1  <=  k  <=  knozzle  and  1  <=  j  <=  jnozzle) 


Region  2  :  wedge  (if  it  exists  at  the  streamwise  location)  ("k  <  1") 


•  Region  3  :  sidewall  (if  it  exists  at  the  streamwise  location)  (k  <  kshen  and 

inozzle  <  i  <  Ishell  ) 


Region  4  :  Flap  (if  it  exists  at  the  streamwise  location)  (j  <  jnozz|e  and 


^nozzle  <  ^  <  ^shell 


•  Regions  5,  6  and  7  :  external  flow  (k  >  kshel,  or  j  >  jshel|) 

The  indices  in  quotes  indicate  areas  where  no  grid  lines  exist.  Four  distinct  types  of 
inconstant  planes  may  be  identified  This  distinction  is  based  on  the  streamwise  distance 
which  determines  the  existence  of  the  flap,  wedge  or  sidewall  at  each  section.  The  range 
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Figure  3-7:  Drff g rent  Types  of  Cross  Sections 


of  x  values  for  each  type  of  cross  section  is  shown  in  Fig  3-3  (b) 

•  Type  I  cross-section:  This  type  of  cross-section  occurs  between  the  upstream 
boundary  and  the  leading  edge  of  the  wedge  (see  Fig.  3-7  (a)).  Note  that  the 
flap  trailing  edge  is  downstream  of  the  wedge  leading  edge  In  the  physical 
plane  of  Fig.  3-6  (a),  the  line  IJ  merges  with  AB  (region  2  vanishes).  The  line 
GK  represents  the  outer  surface  of  the  flap  and  HM  represents  the  inner 
surface  of  the  flap.  The  lines  CL  and  BK  represent  the  outer  and  inner 
sidewall  faces  respectively.  The  boundary  FED  is  the  freestream  boundary. 

•  Type  II  cross  section:  This  type  of  cross  section  occurs  between  the  leading 
edge  of  the  flap  and  the  trailing  edge  of  the  flap  (see  Fig.  3-7  (b))  The  major 
difference  between  this  type  of  cross  section  and  that  of  Type  I  cross  section 
is  that  region  I  is  finite  and  corresponds  to  the  wedge  cross  section. 

•  Type  III  cross  section:  Beyond  the  flap  trailing  edge  but  before  the  wedge 
trailing  edge,  the  region  (GKHM  in  Fig.  3-6)  is  no  longer  solid  (Fig.  3-7(c)). 
The  region  KLCB  continues  to  represent  the  sidewall. 

•  Type  IV  cross  section:  At  a  section  downstream  of  the  wedge  trailing  edge 
(which  occurs  at  the  same  streamwise  location  as  the  sidewall  trailing  edge  in 
the  experiment),  there  are  no  solid  boundaries  in  the  £  =  constant  planes.  The 
computations  extend  some  distance  (a  few  wedge  boundary  layer 
thicknesses)  beyond  the  wedge  trailing  edge  (subsection  3.7)  to  prevent 
corruption  of  data  in  the  region  of  interest  due  to  the  downstream  boundary 
conditions  outlined  in  subsection  3.8.  Since  there  is  no  experimental  data 
available  downstream  in  the  wake  region,  in  the  computation  we  assumed  for 
simplicity  that  a  flat  plate  abuts  the  wedge  trailing  edge  and  that  the  sidewall 
extends  to  the  computational  downstream  boundary.  Under  this  assumption,  a 
Type  IV  cross  section  is  topologically  similar  to  a  Type  III  section  and  hence, 
the  sidewall  is  shown  dashed  in  Fig.  3-7  (d)  indicating  that  it  exists  only  in  the 
computation. 
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3.4.  Theoretical  Model 


A  detailed  description  of  the  theoretical  model  and  numerical  algorithm  may  be  found  in 
[47]  and  is  omitted  here  for  the  purposes  of  brevity.  The  governing  equations  are  the  3-D 
full  mean  compressible  Navier-Stokes  equations  using  mass  averaged  variables  [48], 
Written  in  strong  conservation  form  [49]  and  employing  the  transformed  variables  ^(x,y,z), 
t|(x,y,z)  and  £(x,y,z),  the  equations  are 


where, 


dU  3  F  3  G  3  H  _ 
—  +  —  +  —  +  —  =  0 
3/  3£  3ti  3C 


U  = 


pv  > 
P  w 

_p« 


(7) 


(8) 


J  =  (9) 

3(x,y,z) 

(u,v,w)  are  the  Cartesian  velocity  components  in  the  (x.y.z)  coordinates.  The  density  p, 
static  pressure  p  and  static  temperature  T  are  related  through  the  equation  of  state  p  = 
pRT  where  R  is  the  gas  constant  and  J  is  the  Jacobian  of  the  transformation.  The  total 
energy  per  unit  mass  is  given  by  e  =  e,  +  0.5(u2+v2+w2),  where  the  internal  energy  per  unit 
mass  e<  is  equal  to  CVT  and  Cv  is  the  specific  heat  at  constant  volume.  The  flux  vectors  F, 
G,  and  H  include  the  effects  of  turbulence  in  the  eddy  viscosity  (e)  formulation.  The 
molecular  dynamic  viscosity  p  is  given  by  Sutherland  s  law.  The  molecular  and  turbulent 
Prandtl  numbers,  Pr  and  Prt,  are  0.73  (air)  and  0.9  respectively. 


In  the  vicinity  of  a  surface,  a  special  subset  of  the  Navier-Stokes  equations  is 
employed  [47],  Two  types  of  sublayer  regions  are  defined:  (1)  "line  sublayer"  and  (2) 
"corner  sublayer"  .  These  two  different  regions  are  distinguished  by  their  relative  size  in 
respect  to  the  principal  radii  of  curvature  of  the  surface.  The  line  sublayer  corresponds  to  a 
sublayer  whose  height  (i.e.,  approximately  the  height  of  the  transition  wall  region  of  the 
turbulent  boundary  layer)  is  small  compared  to  the  principal  radii  of  curvature  of  the 
surface  For  this  category,  the  local  surface  appears  "flat"  within  the  scale  of  the  sublayer 
(i.e.,  surface  curvature  effects  may  be  neglected).  The  "corner  sublayer"  corresponds  to  a 


sublayer  where  the  height  of  the  sublayer  is  comparable  or  larger  than  one  of  the  local 
principal  radii  of  curvature  of  the  surface.  A  sharp  corner  (i.e.,  the  local  intersection  of  two 
planes)  is  a  corner  sublayer  since  one  principal  radius  of  curvature  is  zero.  The  code 
requires  that  the  comers  be  90  degrees,  a  condition  satisfied  by  the  nozzle  under 
consideration. 


For  a  line  sublayer,  the  assumption  is  that  the  effects  of  diffusion  of  momentum  and 
energy  in  the  direction  normal  to  the  surface  are  substantially  greater  than  those  in  the 
plane.  The  governing  equations  reduce  to  a  set  of  three  ordinary  differential 
equations  [44,  47]  which  are  solved  with  the  Keller  Box  Scheme  [50]. 

In  the  vicinity  of  a  corner,  the  comer  sublayer  equations  [44,  47]  are  solved.  These 
equations  are  also  obtained  asymptotically  from  the  mean  compressible  Navier-Stokes 
equations  and  reduce  to  two  coupled  nonlinear  partial  differential  equations.  The  solution 
of  these  equations  is  obtained  with  second-order  accurate  finite  differences  and  solved 
using  Newton's  method. 

The  two-layer  turbulent  eddy  viscosity  model  of  Baldwin  and  Lomax  [51]  is  employed  in 
regions  where  no  mixing  layer  exists.  In  this  model,  the  inner  eddy  viscosity  is  given  by 

e(  =  P0c/O)2o)  (10) 

where  k  =  0.40  is  von  Karman's  constant,  /  is  the  Buleev  length  scale  [52,  53,  54],  D  is  the 
Van  Driest  damping  factor  and  co  is  the  absolute  value  of  the  vorticity  (co=|  V  x  v|) 


The  outer  eddy  viscosity  is  given  by 
£0=P*  C  cpF  wakeF  Kleb 


(ID 


where  k  (  =  0  0168)  and  Ccp  (see  below)  are  constants.  The  outer  function  Fwake  is  given 


by 

and, 


wake 


I  F 

maic  n 


=  'outer  =  'WO 


(12) 


(13) 


and  the  maximum  is  evaluated  over  the  entire  boundary  layer.  The  quantity  /  is  the 


value  or  /  where  the  maximum  value  of  Fouter  occurs.  The  value  of  Klebanoff  mtermittency 
correction  in  Eqn.  1 1  is  given  by 

^xith  ~  (1  +  5.5  (^)6)-'  (14) 


where  Ckte5  is  a  constant. 


Recent  calculations  [55]  have  established  that  the  constants  Ccp  and  Ckleb  are  functions 
of  the  Mach  number  and  the  wake  parameters.  For  the  freestream  conditions  established 
in  Table  3-1,  the  values  of  these  two  constants  were  determined  to  be  1.3  and  0.3, 
respectively. 

The  inner  and  outer  eddy  visosity  models  are  combined  to  form  the  turbulent  eddy 
viscosity  e  in  the  following  manner.  First,  profiles  of  e,  and  e0  are  obtained  on  each 
coordinate  line  emanating  from  the  boundary.  This  implies  that  the  coordinate  lines  are 
approximately  orthogonal  within  the  boundary  layer.  The  first  point  nearest  the  boundary  at 
which  exceeds  eq  is  denoted  the  "cross-over  point"  .  The  turbulent  eddy  viscosity  e  is 
then  equal  to  Ej  for  all  points  between  the  boundary  and  the  cross-over  point,  and  is  equal 
to  Eq  for  all  points  above  and  including  the  cross-over  point.  Further  details  of  the 
implementation  of  the  Baldwin  Lomax  eddy  viscosity  model  in  the  vicinity  of  corners 
(internal  and  external)  are  outlined  in  subsection  3.6. 

For  the  free  shear  layer  region  (see  Fig.  3-4),  the  mixing  layer  model  [11, 56]  is 
employed.  In  this  case, 

E  =  p/2  I  co|  (15) 


where, 


i  =  cmixu/\*>L* 


^  ~~  I  v  l/najc  I  V  I min 


C,  =  0.13 

mix 


The  search  for  Icol^  is  conducted  in  a  region  on  either  side  of  the  location  of  |v| 
(subsection  3.6). 


It  may  be  noted  that  this  set  of  governing  equations  in  conjunction  with  the  numerical 
algorithm  employed  (see  subsection  3.5)  is  extendable  to  unsteady  flows  under  the 
assumption  that  the  characteristic  turbulent  diffusion  time  across  the  computational 


sublayer  i.e., 
‘*ff= 


(18) 


is  small  compared  to  the  characteristic  time  scale  of  the  unsteady  mean  motion  tmean.  In 
the  above  expression,  zm  is  the  height  of  the  computational  sublayer,  which  is  typically  less 
than  or  equal  to  60  wall  units  (i.e.,  zmu./vw  <  60,  where  u.  is  the  local  friction  velocity  and 
vw  is  the  kinematic  viscosity  at  the  wall  [21]).  The  quantities  pr,  er  and  qr  are  characteristic 
values  for  the  density,  molecular  dynamic  viscosity  and  turbulent  eddy  viscosity  in  the 
computational  sublayer  (evaluated,  for  example,  at  the  edge  of  the  sublayer).  By  virtue  of 
the  use  of  the  unsteady  Reynolds-averaged  equations,  we  assumed  that  the  time  scale  of 
the  unsteady  mean  motion  tmean  is  similar  to  the  characteristic  development  time  of  the 
mean  flow  tc,  i.e., 


where 


L 


u 


CO 


where  L  is  a  characteristic  length  scale  of  the  mean  flow  (e  g.,  the  airfoil  chord  in  the  case 
of  the  transonic  airfoil  computations  of  Levy  [35]),  and  u^  is  the  characteristic  velocity  of 
the  mean  flow.  The  use  of  the  computational  sublayer  technique,  therefore,  implies 

‘diff  <<:  ‘c  0r  ‘diff  <<  ‘mean 

Due  to  the  exceedingly  small  height  of  the  computational  sublayer  (typically  one  percent  of 
the  boundary  layer  thickness),  this  constraint  is  generally  met  in  practice. 


3.5.  Numerical  Algorithm 

The  numerical  algorithm  is  a  hybrid  explicit-implicit  technique  which  combines  the 
explicit  finite-difference  algorithm  of  MacCormack  [43,  57]  for  the  full  mean  compressible 
Navier-Stokes  equations  and  the  implicit  finite-difference  Box  Scheme  of  Keller  [50]  for  the 
asymptotic  form  of  the  mean  compressible  Navier-Stokes  equations  in  the  viscous 
sublayer  and  transition  wall  regions  of  the  turbulent  boundary  layers. 

For  the  purposes  of  implementation  of  the  explicit  and  implicit  algorithms,  the 
transformed  domain  is  divided  into  two  types  of  regions  -  ordinary  and  sublayer.  These 
regions  are  shown  in  Fig.  3-8  for  each  type  of  section.  The  ordinary  region  consists  of  grid 
points  updated  by  the  explicit  algorithm  while  the  sublayer  regions  are  resolved  with  the 
implicit  algorithm.  The  grid  points  in  the  sublayer  region  are  generated  separately  from  the 
ordinary  grid  with  geometric  stretching.  The  ordinary  and  sublayer  regions  are  interfaced 
along  the  "line  of  matching  points",  which  is  the  row  of  ordinary  grid  points  adjacent  to  each 
solid  boundary  [47], 

The  explicit  algorithm  of  MacCormack  is  implemented  using  a  symmetric  sequence  of 
time  split  operators.  The  method  is  second-order  accurate  in  both  space  and  time  and  is 
relatively  easy  to  program  in  a  numerical  code.  Further,  it  is  easily  and  efficiently  vectorized 
and  is  robust.  Its  overall  stability  is  relatively  insensitive  to  the  initial  conditions  employed. 
Thus,  the  initial  condition  (typically  a  guess)  need  not  be  particularly  close  to  the  final 
(unknown)  solution  [44],  A  fourth-order  pressure  damping  term  [57]  is  incorporated  to 
prevent  numerical  instability  in  the  presence  of  strong  shock  waves. 

The  implicit  algorithm  is  employed  in  the  viscous  sublayer  and  transition  wall  regions  of 
the  turbulent  boundary  layers.  This  region  is  denoted  the  "computational  sublayer"  or 
"sublayer"  for  brevity.  This  hybrid  algorithm  has  been  proven  to  be  very  efficient.  For  a  3-D 
oblique  shock-turbulent  boundary  layer  interaction  [47],  the  hybrid  algorithm  is  a  factor  of 
16  to  21  times  faster  than  a  vectorized  version  of  MacCormack's  explicit  algorithm  alone 
using  time-split  operators.  The  matching  of  the  explicit  and  implicit  algorithms  is  achieved 
through  the  stress  and  heat  flux  tensors  [47]  and  is  omitted  here  for  the  purposes  of 
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Figure  3-8:  Ordinary  and  Sublayer  Regions 


brevity. 


The  governing  equations  and  the  Baldwin-Lomax  turbulent  eddy  viscosity  model  have 
been  successfully  applied  in  conjunction  with  the  hybrid  explicit  algorithm  above  for  a 
variety  of  applications  including  a  3-D  simulated  inlet  [58],  a  3-D  sharp  fin  [59]  and  a  3-D 
swept  compression  corner  [60]  at  supersonic  Mach  numbers. 


3.6.  Details  of  implementation  of  Eddy  Viscosity 

As  mentioned  in  subsection  3.4,  the  two  layer  eddy  viscosity  model  of  Baldwin  Lomax 
model  is  employed  in  regions  at  close  proximity  to  the  wall,  and  the  mixing  layer  model  is 
employed  in  the  wake  region  beyond  the  flap  trailing  edge.  In  general,  the  Buleev  length 
scale  is  taken  to  be  the  distance  normal  to  the  surface  and  e  is  defined  relative  to  the  local 
surface  boundary  layer.  In  the  vicinity  of  an  internal  corner,  the  eddy  viscosity  (Baldwin 
Lomax)  is  implemented  by  dividing  the  tj-C  plane  into  two  areas  in  a  manner  similar  to  the 
method  of  Hung  and  MacCormack  [61]  and  detailed  in  [47],  The  approach  is  illustrated  in 
Fig.  3-9  for  the  case  of  the  corner  formed  by  the  flap  and  the  sidewall  in  a  section  of  Type 
II.  The  line  dividing  the  two  areas  is  taken  to  be  (for  example) 

C  =  ^22/<-^(t)-W).  (20) 

where 


KL- 1 


Ati  = 


JL-l 


and  JL  and  KL  are  the  total  number  of  points  in  the  t)  and  C  direction  respectively. 

For  all  points  in  Area  I,  the  inner  and  outer  eddy  viscosity  profiles  are  obtained  along  i;  = 
constant  lines,  and  the  eddy  viscosity  is  switched  from  the  inner  to  the  outer  formulation  at 
the  location  where  E;  >  e0.  For  points  in  Area  II,  the  inner  and  outer  eddy  viscosity  profiles 
are  obtained  along  q  =  constant  lines  and  a  similar  procedure  is  employed.  The  dividing 
line  formula  20  is  chosen  for  convenience  of  numerical  implementation.  The  image  of  eqn 
20  in  the  physical  domain  is,  in  general,  a  curved  line.  The  exact  appropriate  shape  of  the 
dividing  line  in  the  physical  domain  is  not  known  in  general  and  care  must  be  taken  to 
ensure  that  a  pathological  shape  does  not  occur. 
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Figure  3-9:  Implementation  of  Eddy  Viscosity  Model  Near  Internal  Corner 
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For  external  corners  that  occur  along  the  line  of  intersection  of  the  external  flap  surface 
and  the  external  sidewall  face  (see  Region  7  in  Fig.  3-6),  the  implementation  of  the 
turbulence  model  is  not  clear  since  the  specification  of  the  value  of  the  Buleev  length  scale 
is  difficult.  In  this  research,  a  simple  engineering  approach  is  taken.  The  model  for  the  eddy 
viscosity  in  Region  7  is  obtained  by  extending  on  the  same  eta  or  zeta  line  to  the  line 
Wsheii+k‘ksheii  (see  Fi9-  3-1  °)-  This  eliminates  the  need  for  defining  a  Buleev  length  scale 
in  this  region. 

In  sections  of  Types  I  and  II,  the  Baldwin-Lomax  model  is  exclusively  employed.  In 
sections  of  Type  III  however,  both  the  Baldwin-Lomax  and  the  mixing  layer  models 
(subsection  3.4)  are  applied  in  different  regions  of  the  flow.  The  areas  demarcating  the 
application  regions  of  the  two  models  are  shown  schematically  in  Fig.  3-11.  In  regions  a,b 
and  e,  the  Baldwin-Lomax  model  is  utilized  (along  j  =  constant  lines,  k  =  constant  lines  (j  < 
inozzie)  and  k  =  constant  lines  (j  >  jshe,|)  respectively.  The  height  hw  (typically  about  one 
8wedge)  is  a  function  of  the  streamwise  distance  and  includes  about  10  ordinary  grid  points. 
In  regions  c  and  d,  the  mixing  layer  model  is  employed.  The  search  for  |w|max  is  conducted 
from  top  of  region  "a"  in  Fig.  3-1 1  to  the  top  of  the  grid  (free  stream  boundary). 

3.7.  Coordinate  Transformation  and  Geometrical  Simplifications 

The  grid  is  generated  in  axis-normal  planes  and  stacked  together.  A  limited  number  of 
geometrical  simplifications  of  a  non-critical  nature  are  made  to  simplify  the  implementation 
of  the  numerical  algorithm  and  boundary  conditions  For  example,  a  blunt  wedge  leading 
edge  (as  in  the  experimental  configuration),  results  in  a  90°  kink  in  the  £=0  lines  in  the 
internal  flow  region  (see  Fig.  3-12)  and  leads  to  problems  in  the  computation  of  E, 
derivatives.  The  wedge  is  therefore  assumed  to  be  sharp  by  running  the  tangent  at  (xw  ,, 
yw  ,)  to  the  centerline  of  the  nozzle  Other  modifications  are  explained  when  necessary. 

The  grid  spacing  is  determined  to  satisfy  the  following  general  constraints  [47] 

•  Near  a  solid  boundary,  the  first  ordinary  grid  point  away  from  the  boundary 
must  satisfy  the  condition  that  z'm*  <  60.  where  zm*  is  the  height  of  the 
sublayer  in  wall  units  as  defined  as  in  subsection  3  4. 
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Near  a  solid  boundary,  the  first  sublayer  point  away  from  the  boundary  must 
satisfy  the  condition  z'+  <  3 


•  The  grid  must  possess  about  15  points  inside  the  boundary  layer  (including 
sublayer  points)  for  adequate  resolution  of  gradients. 

The  grid  is  stretched  in  a  direction  normal  to  the  boundary  to  satisfy  the  above 
constraints  without  employing  a  prohibitively  large  number  of  points.  The  maximum  stretch 
factor  anywhere  in  the  grid  is  limited  to  approximately  1.3.  For  grid  generation  purposes, 
approximate  values  of  the  boundary  layer  thickness  (See  Table  3-2),  friction  velocity  (u.) 
and  vw  are  determined  either  with  semi-empirical  relations  or  with  the  2-D  compressible 
turbulent  boundary  layer  code  of  York  [55].  In  Table  3-2,  the  value  of  83  is  estimated 
whereas  the  others  are  specified.  For  the  computation  8  points  are  employed  in  the 
direction  normal  to  the  surface  in  each  sublayer.  A  brief  description  of  the  determination  of 
grid  spacing  in  each  of  the  three  principal  directions  follows. 

Table  3-2:  Boundary  Layer  Thicknesses 


B.L.# 

Notation 

Thickness  (cm) 

1  (Flap  inner) 

^  fylap  mnerl 

0  62  at  x=140.16  cm 

2  (Sidewall  inner) 

§2  or  8Sldewa,|  jnner 

0  62  at  x=140.16  cm 

3  (Wedge) 

^3  or  ^w:  tlgt 

0.12  at  x=167.00  cm 

4  (Flap  outer) 

^4  or  fylap  outer 

1.69  at  x=1 40.00  cm 

5  (Sidewall  outer) 

^5  ^sidewall  outer 

1 .69  at  x=1 40.00  cm 

1.  X  (streamwise)  grid  spacing:  The  computational  upstream  boundary  is  located  at 
x=140.16  cm  from  the  leading  edge  (see  Fig.  3-2).  The  boundary  layer  thickness  on  the 
sidewall  inner  surface  (52)  is  utilized  to  scale  the  streamwise  grid  spacing.  The  axial 
planes  are  arranged  to  provide  at  least  five  planes  inside  the  pressure  rise  (internal  flow) 
with  increased  concentration  of  planes  in  the  shock-boundary-layer  interaction  region 
immediately  downstream  of  the  flap  trailing  edge  The  last  station  is  placed  about  four  52 
downstream  of  the  wedge  trailing  edge  Depending  upon  the  total  number  of  planes 
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employed,  two  cases  are  identified:  Case  I  with  32  planes  and  Case  II  (proposed)  with  64 
planes.  In  this  report  only  Case  I  is  discussed.  The  x-locations  of  each  plane  are  shown  in 
Fig.  3-13  and  values  are  provided  in  Table  3-3. 

2.  Y  (vertical)  grid  spacing:  The  y  direction  extends  normal  to  the  horizontal  plane  of 
symmetry.  Fig.  3-14  shows  typical  grids  at  three  sections  of  Types  I,  II  and  III  respectively. 
The  total  number  of  points  in  this  direction  (for  both  cases)  is  50  (=  KL),  with  30  points 
inside  the  nozzle  (k^^  =  30)  and  20  points  outside  the  nozzle.  To  avoid  compression  of  £ 
grid  lines  originating  inside  the  nozzle,  the  sidewalls  are  modified  to  be  straight  (horizontal) 
from  the  trailing  edge  of  the  flap.  This  corresponds  to  changing  the  value  of  ys  from  0.606 
cm  to  3.775  cm.  Since  the  specification  of  boundary  conditions  on  the  downstream  edge  of 
the  flap  is  not  clear,  by  construction,  no  grid  lines  originate  there  (this  follows  from  equation 
6).  Consider  the  trailing  edge  of  the  flap  shown  magnified  in  Fig.  3-15).  Based  on  the 
criteria  for  the  first  point  away  from  the  boundary  (see  beginning  of  this  subsection),  in 
general,  Au  ^  flap  thickness  *  A,.  The  y  grid  spacing  is  thus  undesirably  discontinuous  at 


this  section.  To  achieve  continuity  of  the  y-grid  spacing  at  the  flap  trailing  edge,  the  height 
of  the  upper  surface  (yb)  is  modified  from  yb  =  3.856  to  yb  =  3.775  cm.  This  corresponds  to 
choosing  Au'  =  flap  thickness  =  A,’  =  minimum(Au,  A,)  where  the  primed  quantities  are  the 
finally  employed  values.  This  approach  leads  to  a  reduction  of  0.081  cm  in  the  flap  trailing 
edge  thickness  which  is  approximately  5  percent  of  the  external  boundary  layer  thickness. 
The  value  of  p  (see  Fig.  3-3  (c))  changes  from  18.65°  to  19.16°.  For  a  freestream  Mach 
number,  M  =  1.2,  a  simple  Prandtl  Meyer  expansion  analysis  indicates  an  acceptable 

oo 

variation  in  pressure  (2.7%)  and  Mach  number  (1.0%).  This  is  comparable  to  (or  less  than) 
the  uncertainty  introduced  by  the  boundary  layer  displacement  effect  (since  the  initial 
boundary  layer  thickness  is  not  known).  The  height  H  (=  6.8  cm)  representing  the  height  of 
the  external  flow  above  the  flap  at  the  upstream  boundary  is  chosen  to  ensure  that  the 


Mach  waves  arising  from  the  flap,  after  reflection  from  the  outer  boundary  (due  to  the 
freestream  boundary  conditions  applied),  pass  through  the  downstream  boundary  and  do 
not  pollute  the  computation  in  the  region  of  interest. 
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a)  Actual  (Experiment) 


b)  Modified  (Computation) 

Figure  3-12:  Modification  of  Wedge  Leading  Edge 


In  the  internal  flow  region  (j  <  jnozzle,  k  <  kn0ZZ|e),  the  y  grid  is  stretched  to  concentrate 
more  points  near  the  wedge  and  flap  inner  surfaces.  In  region  5,  the  y  spacing  is  adjusted 
to  concentrate  grid  points  only  at  the  flap  outer  surface.  The  y  spacing  in  region  5  is 
extended  to  region  7.  In  region  6  (j  >  jshel)l  k  <  kshen),  the  grid  is  stretched  away  from  k  » 
kshel,  towards  k  =  1  ending  in  several  equal  spaced  grid  lines.  The  grid  spacing  at  the 
junction  of  regions  6  and  7  varies  in  a  smooth  fashion. 

3)  Z  (spanwise)  grid  spacing  Fig  3-13  shows  typical  j  =  constant  planes  of  the  grid  A 
total  of  39  points  are  employed  in  the  z  direction  (JL  =  39),  with  21  (=  Jnozzie)  points  in  the 
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Table  3-3:  X-locations  (cm.) 


Stn. 

X  (cm. ) 

Stn . 

X  (cm. ) 

Stn . 

X (cm. ) 

Stn . 

X  (cm. ) 

1 

140.1600 

9 

150.8829 

17 

157.2641 

25 

162.0490 

2 

141.7600 

10 

151.8758 

18 

157.8506 

26 

162.8323 

3 

143.3341 

11 

152.8056 

19 

158.3998 

27 

163.6783 

4 

144.8079 

12 

153.6761 

20 

158.9224 

28 

164.5920 

5 

146.1879 

13 

154.4911 

21 

159.4552 

29 

165.5790 

6 

147.4799 

14 

155.2542 

22 

160.0308 

30 

166.6244 

7 

148.6897 

15 

155.9688 

23 

160.6524 

31 

167.6698 

8 

149.8223 

16 

156.6377 

24 

161.3238 

32 

168.7153 

interior  nozzle  flow.  To  obtain  continuity  of  z-grid  spacing  adjacent  to  the  wall,  for  reasons 
similar  to  those  at  the  trailing  edge  of  the  flap,  the  sidewall  exterior  surface  is  assumed 
parallel  to  the  x-axis  i.e.,  the  4.75°  inward  slant  is  ignored.  A  Prandtl  Meyer  analysis  similar 
to  that  carried  out  for  the  flap  thickness  indicates  that  the  effects  are  modest.  Uniformity  is 
achieved  by  assuming  that  the  sidewall  thickness  is  arbitrarily  small  and  equal  to  the  value 
necessary  for  adequate  resolution  of  the  interior  and  exterior  side  boundary  layers 

3.8.  Initial  and  Boundary  Conditions 

The  freestream  and  wall  (temperature)  conditions  (from  experiment)  are  summarized  in 
Table  3-1 .  For  the  internal  flow,  the  choice  of  nozzle  pressure  ratio  (NPR  =  2)  is  dictated  by 
two  factors: 

•  throat  is  choked  (NPR  greater  than  1.89),  and 

•  shock  must  hit  wedge  as  far  upstream  as  possible  to  avoid  separation  from 
extending  into  the  wake  region  where  the  mixing  layer  eddy  viscosity  model  is 
employed. 

The  boundary  conditions  employed  may  be  classified  into  the  following  types: 

•  Upstream  boundary 

•  Downstream  boundary 


she! 


•  Solid  boundaries 

•  Symmetry  boundaries 

•  Freestream  Boundaries 


Upstream  (Inflow)  Boundary:  The  upstream  boundary  profile  may  be  classified  into  two 
regions:  external  flow  and  internal  flow.  In  the  external  flow  regions,  the  flow  profiles  in 
regions  5  and  6  are  specified  by  assuming  that  the  boundary  layers  (numbers  4  and  5 
respectively)  develop  from  the  tip  of  the  nozzle  (Fig.  3-2).  Under  this  assumption,  the  flow 
profile  is  generated  separately  with  a  2-D  compressible  boundary  layer  code  [55).  In  region 
7,  the  boundary  layer  profiles  are  combined  in  a  manner  similar  to  that  employed  for  the 
eddy  viscosity  i.e.,  points  below  the  dividing  line  in  region  7  obtain  the  flow  variable  values 
from  region  6  along  j  =  constant  lines  while  points  above  the  line  obtain  their  values  from 
region  5  along  k  =  constant  lines. 


In  the  internal  flow  region  the  specification  of  the  boundary  condition  is  somewhat  more 
complex.  The  application  may  be  divided  into  inviscid  and  viscous  regions.  In  the  inviscid 


region  (where  the  flow  is  subsonic),  from  the  theory  of  characteristics  [62], 
dt  dx 


The  values  of  pressure  at  the  previous  iteration  permit  the  computation  of  the  value  of 
pressure  at  the  next  time-step  using  one-sided  first-order  accurate  differences  in  space 
and  time.  Since  the  total  mass  flow  rate  depends  upon  the  total  pressure,  total  temperature 
and  the  throat  area,  (Wj  ~  P/V7j),  the  total  pressure  and  total  temperature  are  held  fixed.  In 
the  inviscid  region,  the  total  pressure  and  the  pressure  (at  timestep  (n+1))  yield  the  Mach 
number  (and  hence  the  other  flow  parameters)  at  (n+1).  In  the  viscous  region,  the  velocity 
profile  in  the  boundary  layer  may  be  obtained  from  the  incompressible  law  of  the 


wall/wake  [63,  55),  with  M  «  1 


,  .  1 ,  /u»,  n  7U  i,ny, , 

!  =  u{-ln(- — )  +  B  +  — .w/r(-£) 


(22) 
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In  equation  22,  u.  denotes  the  friction  velocity  (known  since  the  magnitude  of  the  wall 
shear  stress  is  maintained  fixed),  B  is  a  constant  (=5  1),  n  is  the  wake  strength  parameter 
(see  equation  24),  The  magnitude  of  the  wall  shear  stress  tw  (c,  =  6  1  X  10  3  based  on  the 
nominal  Mach  number),  the  boundary  layer  thickness  (5)  and  total  temperature  (Tt)  are 
maintained  fixed.  In  the  above  formula  therefore, 

(23> 


>(—)-/?) 
k  v 


Equations  22  through  24  permit  the  determination  of  velocity  u  in  the  boundary  layer. 

The  temperature  may  then  be  computed  with  the  formula  : 

T=T,-J-u 2  ;  v  =  iv  =  0  (25) 

ZLp 

The  upstream  boundary  condition  in  the  internal  flow  region  therefore  varies  with  the 
computation.  We  found  this  boundary  condition  quite  satisfactory  as  is  evident  from  results 
presented  in  subsection  3.1 1 . 

Downstream  (Outflow)  Boundary:  The  outgoing  flow  is  supersonic,  and  the  outflow 
boundary  is  essentially  perpendicular  to  the  outgoing  flow.  The  flow  characteristics 
therefore  point  from  the  inside  of  the  computational  domain  to  the  outside,  and  thus  the 
outflow  boundary  flow  variables  must  be  determined  from  the  interior  flow  solution.  The 
conventional  zero-gradient  extrapolation  d/d£  =  0  condition  is  employed.  To  prevent  the 
effects  of  this  boundary  condition  from  affecting  the  region  of  interest,  the  computational 
domain  is  extended  two  ^-planes  downstream  of  the  trailing  edge  of  the  wedge. 

Solid  boundaries:  At  a  solid  boundary,  the  following  boundary  conditions  are  applied 
p  v  =  0 

^  =  0  (26) 


where  v=  (u,v,w)  is  the  Cartesian  velocity  vector,  p  is  the  static  pressure,  the  direction  n  is 
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normal  to  the  boundary  and  Tw  is  a  specified  wall  temperature  (see  Table  3-1)  The 
boundary  condition  for  static  pressure  is  an  approximation  to  the  exact  boundary  condition 
derivable  from  the  normal  component  of  the  momentum  equations,  and  has  been 
succesfully  employed  for  a  variety  of  flows  exhibiting  shock-boundary  layer 
interaction  [15,  22,  44,  64,  65] 

Symmetry  boundaries.  On  a  symmetry  boundary  the  following  conditions  apply  : 

-  A 

p  v .n  =  0 

A(pvxn)  =  0 

on 

(27) 

^  =  0 

Bn 

Z.o 

dn 


where  n  is  the  direction  normal  to  the  symmetry  boundary. 


Freestream  boundaries:  At  these  boundaries,  the  flow  variables  are  obtained  from  the 

conditions  outlined  in  Table  3-1  and  are  applied  in  the  form  : 
p  =  specified 

u  =  specified 

(28) 

v  =  w  =  0 
pe  =  specified 

Fig.  3-16  shows  schematically  the  boundary  conditions  as  applied  at  each  of  the  three 
types  of  cross  sections.  In  all  computations  described  in  this  research,  the  flow  initial 
condition  is  obtained  by  simply  propagating  the  profiles  at  the  upstream  boundary  (c  0)  to 
all  downstream  grid  points. 


3.9.  Computational  Details 

The  marching  of  the  solution  toward  steady-state  is  carried  out  in  many  phases.  For  the 
first  approximately  0.1  Tcext  (characteristic  time  based  on  external  flow  -  the  time  required 
for  a  fluid  particle  to  travel  from  the  upstream  to  the  downstream  end  of  the  physical 
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Figure  3-16:  Boundary  Conditions  at  Each  Type  ot  Section 
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domain  in  the  external  flow  inviscid  region),  the  sublayer  solution  is  updated  every  explicit 
iteration,  This  is  necessary  because  the  initial  condition  is  obtained  by  propagating  the 
upstream  boundary  condition  to  all  streamwise  planes  and  is  thus  far  removed  from  the 
steady-state  solution.  Therefore,  the  explicit  solution  develops  quite  rapidly  and  the 
sublayer  solution  must  be  updated  every  explicit  step.  In  the  initial  stages  also,  a  "transient 
damping  procedure"  is  employed  to  prevent  numerical  instability  of  the  hybrid  algorithm.  A 
"poor"  initial  guess  often  leads  to  transient  spatial  oscillations  in  static  pressure.  These 
oscillations  disappear  after  approximately  one  characteristic  time.  To  prevent  these 
unphvsical  transient  oscillations  from  adversely  affecting  the  line  sublayer,  the  transverse 
pressure  gradient  3 p/dy'  in  the  line  sublayer  is  set  to  zero  for  an  initial  time  interval 
approximately  equal  to  one  Tcext.  In  the  present  instance,  numerical  instability  is  observed 
(originating  at  the  external  corner)  when  the  transverse  pressure  gradient  is  switched  on 
even  after  one  Tc  ext.  The  problem  is  "solved"  by  computing  this  gradient  normally  at  all 
points  except  near  the  external  corner  (see  Fig.  3-17)  where,  over  the  last  few  grid  points 
the  value  of  the  gradient  is  arbitrarily  set  to  a  value  determined  by  a  linear  profile  to  obtain 
zero  value  at  the  corner.  The  effect  of  this  assumption  is  not  expected  to  be  significant. 

After  about  one  characteristic  time  Tcext,  the  sublayer  solution  develops  at  a  slower 
pace  and,  since  the  sublayer  solution  takes  several  times  the  CPU  as  the  explicit  solution, 
it  is  updated  only  once  about  every  eight  time  steps  and  the  solution  marched  to 
convergence. 

Convergence  is  determined  by  checking  flow  variation  over  roughly  one  Tcext.  Typically, 
the  solution  vector  U  (Eqn.  8)  and  the  eddy  viscosity  e  are  monitored  on  the  ordinary  grid, 
while  the  velocities,  temperatures,  their  gradients  and  e  are  monitored  on  the  sublayer  and 
corner  grids.  Since  certain  quantities  (e  in  the  freestream  for  example)  obtain  very  small 
values,  relative  changes  in  these  may  be  rather  large  and  misleading.  Variations  in  all 
quantities  are  therefore  monitored  relative  to  local  values  as  well  as  over  global  reference 
values  to  determine  convergence.  One  interesting  observation  is  that  the  convergence 
properties  of  the  sublayer  algorithm  improve  significantly  when  the  eddy  viscosity  (in  the 
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Figure  3-17:  Computation  of  Transverse  Pressure  Gradient  in  Sublayer 

at  External  Comer 

sublayer)  is  averaged  over  successive  sublayer  solutions.  The  cause  of  this  is  presently 
unclear.  Convergence  is  achieved  over  most1  parts  of  the  domain  by  marching  the  explicit 
solution  forward  by  approximately  4.10  Tc  ext  flow  development  time  corresponding  to  1694 
iterations  on  the  ordinary  grid.  This  requires  roughly  13  hours  CPU  on  the  2-pipe  CYBER 
205  at  the  John  von  Neumann  Center  (JVNC). 


'see  subsection  3.11  regarding  possible  existence  of  pressure  waves  in  the  vicinity  of  the  sidewall -wedge 
comer 
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3.10.  Comparison  with  Experiment 

The  accuracy  of  the  numerical  results  are  validated  by  comparison  with  experimental 
static-pressure  data  of  Mason  and  Abeyounis  [46],  The  locations  of  the  14  rows  of  static 
pressure  orifices  in  the  experimental  configuration  are  shown  in  Fig.  3-18  and  summarized 
in  Table  3-4.  The  computed  flowfield  is  observed  to  converge  to  a  steady  state  at  all 
streamwise  planes  over  a  large  spanwise  (t})  portion  of  the  domain  from  the  symmetry 
plane  ABFE  (Fig.  3-5)  to  a  spanwise  station  near  the  sidewall.  In  the  vicinity  of  the  corner 
formed  by  the  wedge  and  the  sidewall,  however,  an  unsteady  solution  is  obtained  as 
described  in  subsection  3.11.  The  experimental  data  includes  mean  surface  pressure 
measurements  only.  Consequently,  only  the  static  pressure  distributions  on  the  centerlines 
of  the  flap  (inner  surface),  flap  (outer  surface)  and  wedge  are  compared  in  the  present 
report.  These  correspond  to  rows  1  (Fig.  3-18  (a)),  4  (Fig  3-18  (b))  and  10  (Fig.  3-18  (e)) 
respectively.  In  the  following  graphs,  the  static  pressure  is  non-dimensionalized  with  the 
total  pressure  in  the  nozzle  at  the  upstream  station  and  the  dimension  d,  the  maximum 
height  of  the  model  (Fig.  3-3)  is  employed  to  normalize  the  axial  distance  x. 

The  pressure  along  the  wedge  centerline  is  plotted  in  Fig.  3-19.  As  observed  in 
experiment,  the  flow  expands  as  the  flap  trailing  edge  (x/d  =  1.0)  is  reached.  Downstream 
of  the  flap  trailing  edge  (nozzle  exit),  the  static  pressure  rises  somewhat  abruptly  indicating 
the  possible  presence  of  a  shock  on  the  surface  of  the  wedge.  Beyond  this  pressure  rise 
(x/d  >  1.15),  experimental  observations  indicate  a  roughly  constant  static  pressure  profile. 
From  Fig.  3-19,  the  computations  overpredict  the  static  pressure  downstream  of  the  wedge 
shock  by  approximately  5  percent.  The  conclusion  made  by  the  experimental  investigators 
regarding  the  existence  of  flow  separation  on  the  wedge  (beyond  the  shock)  is  to  be 
confirmed  with  particle  tracing  techniques  (subsection  3.11). 

The  pressure  along  the  flap  internal  centerline  is  plotted  in  Fig.  3-20.  Note  that  (he  last 
experimental  station  is  located  at  roughly  x/d  =  1.0  which  corresponds  to  the  flap  trailing 
edge.  The  pressure  drops  as  the  flap  trailing  edge  is  reached.  The  calculated  and 
experimental  static  pressures  are  in  good  agreement  over  the  entire  length  of  the  flap.  The 


(a)  Internal  flap  pressure  orifices. 


(b)  External  flap  pressure  orifices. 


(c)  External  sidewall  (right  sidewall) 
pressure  orifice*. 


Left 

(d)  Internal  sidewall  (left  sidewall) 
pressure  orifices. 


(e)  Wedge  pressure  orifices. 

Figure  3-18:  Location  of  Experimental  Static  Pressure  Orifices. 

From  [46] 


Table  3-4:  Summary  of  Static  Pressure  Orifice  Locations 


Row  # 

Location 

Figure 

1,2,3 

Internal  Flap 

3-18  (a) 

4,5 

External  Flap 

3-18  (b) 

6,7,8 

External  Sidewall 

3-18  (c) 

9 

Internal  Sidewall 

3-18  (d) 

10-14 

Wedge 

3-18  (e) 

assumption  made  on  the  flap  trailing  edge  thickness  apparently  did  not  affect  the 
computations. 

The  pressure  coefficient  along  the  flap  external  centerline  is  plotted  in  Fig.  3-21.  The 
comparison  shows  reasonable  agreement  for  x/d  <  0.45.  For  values  of  x/d  >  0.45  however, 
the  computed  pressure  coefficient  is  somewhat  higher  than  observed  in  experiment.  The 
cause  of  this  discrepancy  is  under  investigation.  For  the  low-expansion-ratio  configuration, 
the  experimental  static  pressure  shows  spanwise  variation  beyond  the  nozzle  (x/d  >  0  8) 
on  the  flap  outer  surface  [46]  and  further  comparison  is  necessary 

3.11.  Analysis  of  Computed  Flow  Field  Structure 

One  of  the  many  advantages  of  numerical  simulation  of  flow  fields  is  the  capability  of 
complete  analysis.  In  particular,  it  is  perhaps  easier  to  detect  the  existence  of  critical  flow 
structures  leading  to  a  more  complete  understanding  of  the  salient  features  of  flows  where 
experimental  visualization  is  difficult. 

Three-dimensional  computations  as  described  above  necessarily  result  in  sizeable 
datasets.  The  processing  of  such  large  sets  of  numbers  is  a  non-trivial  task  and,  for  the 
current  computation,  is  currently  under  progress.  The  primary  focus  of  the  flow  field 
analysis  is  on  the  study  of  particle  traces  generated  with  the  computed  flow  field.  In  the 
present  instance,  it  is  proposed  to  employ  a  3-D  particle  tracing  code  [66].  A  number  of 
particles  are  released  in  different  parts  of  the  domain  and  their  motion  is  studied  to  detect 
coherent  flow  structures  such  as  vortices  and  regions  of  separation.  It  is  proposed  to 
analyze  the  paths  of  these  particles  with  sophisticated  color  graphics  software 
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Figure  3-19:  Static  Pressure  Variation  on  Wedge  Centerline  - 
Comparison  with  Experiment 
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Figure  3*20:  Static  Pressure  Variation  on  Flap  Inner  Centerline 
Comparison  with  Experiment 


Figure  3-21 :  Static  Pressure  Variation  on  Flap  Outer  Centerline  - 
Comparison  with  Experiment 
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(MOVIE. BYU)  to  provide  a  proper  perspective. 

An  indication  of  the  appropriateness  of  the  upstream  boundary  condition  for  internal  flow 
is  obtained  from  Fig.  3-22  where  the  upstream  static  pressure  is  plotted  against  the 
number  of  iterations.  The  pressure  displays  oscillatory  behavior  as  the  solution  develops. 

The  amplitude  of  these  oscillations  is  quite  small  (-  0.3  %).  These  oscillations  persist  even 
at  "steady  state"  conditions  and  do  not  affect  the  flow  solution  to  any  significant  extent 
insofar  as  the  surface  pressures  are  concerned. 

Preliminary  analysis  of  the  results  indicates  the  existence  of  pressure  waves  at  sidewall- 
wedge  corner  where  an  unsteady  solution  is  observed.  This  conclusion  is  supported  by 
Fig.  3-23  which  shows  the  development  of  pressure  along  an  axial  line  on  the  wedge 
surface  near  the  wedge-sidewall  corner.  The  peaks  marked  A,B  and  C  in  Fig.  3-23 
demonstrate  the  existence  of  a  pressure  wave.  It  is  not  clear  now  if  these  are  actual 
physical  phenomena  or  if  they  are  consequences  of  the  numerical  solution  process  (for 
example  due  to  the  downstream  boundary  conditions).  The  phenomenon  is  being  analyzed 
further. 

Deductions  from  values  of  shear  stresses  also  point  to  the  existence  of  regions  of  ! 

separated  flow  on  the  wedge  (x/d  -  1.0).  Separation  is  also  indicated  on  the  external 

corner  formed  by  the  flap  and  the  sidewall  though  it  is  not  clear  now  if  this  is  due  to  the  ‘ 

assumptions  regarding  the  geometry  of  the  nozzle  sidewall.  j 

I 

3.12.  Conclusion 

Preliminary  results  of  a  complex  nozzle  flow  simulation  indicate  very  good  agreement 

between  theory  and  experiment  along  the  centerline  of  the  nozzle.  Modest  discrepancies  - 

■ 

are  observed  on  the  flap  external  centerline.  The  flow  seems  rather  complex  with  the  ' 

existence  of  a  pressure  wave  system  in  the  wedge-sidewall  corner.  As  indicated  above,  „ 

we  expect  further  analysis  of  the  computed  flow  field  to  result  in  the  detection  of  further  t 

) 

coherent  structures.  Thus  we  should  have  a  better  understanding  of  the  flow  dynamics.  \- 


69 


960  0 


1200.0 


0.7  0.8  0.9  1.0  1.1  1.2  1.3  1.4  1.5  1.8 

X  /  d 


Figure  3-23:  Static  Pressure  Near  Sidewall-Wedge  Corner 
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4.  Conclusion  and  Future  Work 

A  three-dimensional  graphics  based  interactive  method  of  grid  generation  is  developed. 
The  approach  employs  two  phases.  In  the  first  phase,  the  algebraic  multisurface  technique 
of  Eiseman  is  utilized,  together  with  the  precise  point  distribution  technique  of  Smith  et  al. 
In  the  second  phase,  a  series  of  direct  grid  modification  techniques  are  utilized  to  meet 
specific  grid  requirements  not  easily  satisfied  during  the  first  phase.  The  generality  of  the 
method  is  demonstrated  for  a  twin-jet  axisymmetric  nozzle  and  an  aircraft  section. 

A  code  based  on  the  hybrid  explicit-implicit  algorithm  [47]  is  modified  in  this  research  to 
compute  nozzle  flows.  The  principal  modifications  include  the  ability  to  simulate  flow  fields 
for  geometries  with  internal  and  external  flows  and  the  use  of  the  mixing  layer  model  where 
these  flows  meet.  The  modifications  retain  the  vectorization  philosophy  of  the  original  code. 
The  major  coding  modifications  involve  the  solution  of  the  sublayer  equations  on  a  larger 
number  of  surfaces  as  a  result  of  the  more  complex  nozzle  topology  (as  compared  to  the 
3-D  sharp  fin  topology  for  example)  and  the  interfacing  of  the  two  different  eddy  viscosity 
models  employed. 

The  flow  in  two  non-axisymmetric  nozzles  is  computed  under  the  assumptions  of 
horizontal  and  vertical  symmetry.  The  freestream  Mach  number  of  the  external  flow  is  1.2 
and  the  nozzle  pressure  of  the  internal  flow  is  2.0.  One  salient  feature  of  the  computation  is 
the  employment  of  a  different  upstream  boundary  condition  for  the  internal  flow.  This 
condition  depends  upon  the  theory  of  characteristics  tor  inviscid  regions  and  the  law  of  the 
wall/wake  for  the  viscous  regions. 

The  solution  for  the  flow  field  on  the  centerline  is  converged.  Comparison  of  static 
pressure  values  at  the  flap  (inner  and  outer)  and  wedge  surface  centerlines  with  the 
experimental  observations  of  Mason  and  Abeyounis  [46]  indicate  very  good  agreement 
with  modest  discrepancies  in  the  shear  layer  region  where  the  computations  predict 
generally  higher  pressures. 

A  preliminary  analysis  of  the  computed  flow  field  indicates  the  existence  of  a  shock  wave 
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on  the  wedge  in  the  region  beyond  the  flap  trailing  edge.  The  solution  in  the  corner  formed 
by  the  wedge  and  the  sidewall  displays  unsteady  behaviour  (Fig.  3-23).  The  wave  structure 
of  the  flow  in  this  region  may  not  have  been  detected  in  the  experiment  due  to  inadequate 
frequency  response  characteristics  of  the  experimental  apparatus. 

Future  work  will  focus  on  a  more  thorough  analysis  of  the  unsteady  phenomenon 
discussed  above.  It  is  envisaged  that  the  flow  will  be  studied  in  detail  with  numerical 
particle  tracing  techniques  and  color  graphics  to  identify  important  flow  structures  such  as 
areas  of  separation  and  vortices.  A  refined  case  with  more  (64)  axial  planes  is  also  under 
consideration. 
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Nomenclature 


Symbols  relevant  to  the  geometric  configuration  in  the  experiment  are  provided  in  Fig.  3-3. 
Ak  Reciprocal  of  Gk  (Equation  1) 


Baldwin  Lomax  turbulent  eddy  viscosity  model 


Specific  heat  at  constant  pressure 


Specific  heat  at  constant  volume 


Turbulence  model  (B-L)  constant  (Equation  11) 


Turbulence  model  (B-L)  constant  (Equation  14) 


Turbulence  model  (mixing  layer)  constant  (=  0.13)  (Equation  16) 


Pressure  coefficient 


van  Driest  damping  factor  (Equation  10) 


Total  energy  per  unit  mass 


Internal  energy  per  unit  mass 


Flux  vector  in  4  direction  (Equation  7) 


Klebanoff  intermittency  function  (Equation  11) 


Velocity  scale  in  outer  turbulence  model  (Equation  12) 


Outer  function  in  outer  turbulence  model  (Equation  13) 


Wake  function  in  outer  turbulence  model  (Equation  12) 


Flux  vector  in  q  direction  (Equation  7) 


Flux  vector  in  C,  direction  (Equation  7) 


Integrals  of  interpolants  y  (Equation  1 ) 


Index  in  the  %  direction 


Jacobian  of  transformation  (Equation  9) 


Index  in  the  q  direction 


Index  in  the  £  direction;  thermal  conductivity; 


j  (spanwise)  index  corresponding  to  sidewall  inner  surface 


j  index  corresponding  to  sidewall  outer  surface 


Constant  in  outer  eddy  viscosity  formulation 


k  index  corresponding  to  flap  inner  surface 


k  index  corresponding  to  flap  outer  surface 

Buleev  length  scale(B-L)  (Equation  10);  length  scale  (mixing  layer 
Equation  15) 

Freestream  Mach  Number  (Table  3-1) 

Total  number  of  surfaces  in  multisurface  technique  (Equation  1) 
Nozzle  Pressure  Ratio  (Table  3-1) 

Reynolds  number  (Table  3-1) 

Unit  normal  at  surface  (points  inwards) 

General  multisurface  transformation  (Equation  1) 

Parameterized  surfaces  (Equation  1) 

Average  tunnel  total  pressure  (Table  3-1) 

Molecular  Prandtl  number 
Turbulent  Prandtl  number 
Static  pressure 
Gas  constant  for  air 

Spanning  variable  in  multisurface  technique  (Equation  1) 

Skewness  measure  (Equation  5) 

Static  temperature 
Characteristic  time 

Characteristic  time  based  on  external  flow 
Tunnel  total  temperature  (Table  3-1) 

Average  jet  total  temperature  (Table  3-1) 

Wall  temperature  (Table  3-1) 

Time 

Characteristic  velocity  of  mean  flow  (Equation  19) 

Characteristic  turbulent  diffusion  time  across  sublayer  (Equation  18) 
Time  scale  of  unsteady  mean  motion 
Solution  vector  <p,  pu,  pv,  pw,  pe>  (Eqn.  8) 

Velocity  vector  (u,v,w) 

Cartesian  velocity  in  x  direction 


u. 


w 

wi 

X 

y 

z 


z’ 

zm,zm 


GREEK  SYMBOLS 


Friction  velocity  (Equation  23) 

Cartesian  velocity  in  y  direction 
Magnitude  of  maximum  velocity  (mixing  layer) 

(Equation  17) 

Magnitude  of  minimum  velocity  (mixing  layer) (Equation  17) 

Cartesian  velocity  in  z  direction 

Theoretical  mass  flow  rate  in  nozzle 

Axis  in  the  general  streamwise  direction 

Axis  normal  to  the  horizontal  plane  of  symmetry 

Axis  in  the  general  spanwise  direction 

Local  sublayer  distance  normal  to  surface 

Height  of  computational  sublayer 

Sublayer  height  in  wall  units  (non-dimensional) 


4 

\ 

Au 

8 

^  1  ’  lap  inner 
^2’^sidewall  inner 
^3'Awedge 
^4'^flap  outer 

^5'^sidewall  outer 
E 

£i 

eo 

er 

1} 


Maximum  permissible  sublayer  height  on  flap  lower  surface  (Fig.  3-15) 
Maximum  permissible  sublayer  height  on  flap  upper  surface  (Fig.  3-15) 
Actual  sublayer  height  employed  on  flap  lower  surface  (Fig.  3-15) 
Actual  sublayer  height  employed  on  flap  upper  surface  (Fig.  3-15) 
Boundary  layer  thickness 
Inner  flap  5  (Table  3-2) 

Inner  sidewall  8  (Table  3-2) 

Wedge  8  (Table  3-2) 

Outer  flap  8  (Table  3-2) 

Outer  sidewall  8  (Table  3-2) 

Eddy  viscosity 

Inner  eddy  viscosity  (Equation  10) 

Outer  eddy  viscosity  (Equation  1 1 ) 

Characteristic  eddy  viscosity  in  sublayer  (Equation  18) 

Computational  axis  (index  -  j) 
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von  Karman  constant  (  =  0.4)  (Equation  10) 

Molecular  dynamic  viscosity 

Characteristic  molecular  dynamic  viscosity  in  sublayer  (Equation  18) 
Kinematic  viscosity 
Kinematic  viscosity  at  wall 
Wake  parameter  (Equation  24) 

Density 

Characteristic  density  in  sublayer  (Equation  18) 

Angles  (degrees)  surrounding  a  grid  point.  (Equation  5) 

Shear  stress 
Shear  stress  at  wall 
Vorticity 

Magnitude  of  maximum  vorticity  (mixing  layer  model)  (Equation  16) 
Computational  axis  (index  -  i) 

Interpolants  for  multisurface  technique  (Equation  2) 

Computational  axis  (index  -  k) 

jet 

total  (temperature  for  example) 
evaluated  at  wall 
Evaluated  at  freestream 
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